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ing a metallic surface is solved in the wide-band limit. Equations for the time-evolution of the 
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I. INTRODUCTION 



Great strides have been made in recent years in our understanding of the adsorption and reaction of molecules on 
surfaces. Much of this progress is due to the development of electronic structure methods based on density functional 
theory (DFT), from which the potential energy surface (PES) for the molecule-surface interaction can be obtained. 
When combined with classical or quantum dynamical calculations, the ab-initio PES provides a clear picture of the 
molecular pathways and the making and breaking of bonds^. It is well known that DFT-based calculations are 
approximate in their use of a model for the exchange-correlation potential. However, for any dynamical process there 
is another key approximation made. DFT calculations are firmly rooted in the Born-Oppenheimer approximation 
ON ' (BOA); the electronic system remains in the ground state throughout the molecule-surface encounter and all the 
' dynamics takes place on the adiabatic PES 2 . However, the BOA cannot strictly be valid for interactions involving a 
i metallic substrate. This is because there is a continuum of electronic states at the Fermi level, and consequently any 
■ dynamical process will lead to the excitation of electron-hole pairs^. 

How large is the coupling between the nuclear motion and the electronic system? Some very recent studies have 
suggested that non-adiabatic effects are small for closed shell molecules. Dynamical calculations for H2/Pt£ and 
N2/P1U- based on the ground state PES show good agreement with molecular beam scattering experiments, suggesting 
that coupling to electron-holes pairs is not significant for these systems. However, there is strong evidence in other cases 
that non-adiabatic effects are large. Two recent series of experiments have provided striking evidence of electronic 
excitations in surface reactions. First, Nienhaus and co-worker o 5 i 6 i 7 ' 8 adsorbed a range of atomic and molecular 
• • ■ species on thin metal films of silver and copper which made up one contact of a Schottky diode. The hot electrons 



or holes resulting from the dissipation of the chemisorption energy were detected as a chemically induced current, 
or chemicurrent. Second, White and co-worker o 9 ' 10 investigated the scattering of vibrationally excited NO molecules 
from a low work function, caesium-doped gold surface. For molecular beams with a vibrational energy greater than the 
work function the emission of exo-electrons was observed, showing that electronic excitation is a significant channel 
for the dissipation of vibrational energy. 

Most theoretical treatments of non-adiabatic effects use a "nearly-adiabatic" approximation in which the pertur- 
bation of the electronic system is assumed to be weak and slow. This leads to a friction-based description of the 
energy transfer between the adsorbate and the substrate (seeii and references therein) . The friction coefficient can be 
calculated using ab-initio methods^ 2 -, and has been applied to vibrational dampingii and desorption dynamics^. By 
making a connection between the friction description and the forced oscillator model it is also possible to obtain the 
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spectrum of electron-hole pair excitation o 11 ' 15 . However, the friction-based description is not sufficient for strongly 
non-adiabatic effects, and it fails completely for the adsorption of spin-polarised species. Trail et a h 11 ' 15 attempted to 
use the nearly-adiabatic approximation to model the excitation of chemicurrents in the adsorption of H atoms on Cu. 
They found that the friction coefficient diverges at the point in the trajectory where the ground state has a transition 
from being spin polarised (H-atom far from the surface) to non-polarised (H close to the surface). Non-adiabatic 
effects in this case are large and a theoretical description must go beyond the friction approximation. 

In previous wor k 16 ' 17 we have attempted to understand the origin of this divergence in the friction coefficient by 
modelling a strongly non-adiabatic coupling between an adsorbate and a metal surface. We used the simplest possible 
model of a gas-surface interaction: the Newns-Anderson mode l 18 ' 19 , where a single adsorbate level interacts with a 
wide electronic band. Within the mean-field approximation we showed that the time dependence of the occupancy 
of the adsorbate state can be found analytically, and we derived expressions for the rate of energy transfer to the 
surface^ 7 -. The Newns-Anderson model shows the same ground-state spin transition as the ab-initio DFT calculations, 
but the singularity in the energy transfer is removed in the fully non-adiabatic solution. In our previous work we did 
not solve for the evolution of the substrate states, and so we could not derive the spectrum of electronic excitations 
which would be required, for example, in modelling chemicurrents. Our aim in this paper is therefore to extend our 
previous work to a full solution of the non-adiabatic dynamics. In section[H]wc derive an expression for the distribution 
of occupied electronic states and how this evolves with time. The equivalent distribution without electronic excitation 
is considered in section IIII1 which allows us to calculate the spectrum of excitations. Section IIVI gives an outline of 
the methodology we have used to compute the excitation spectrum. We demonstrate the behaviour of the model in 
section [Vj and concluding remarks are made in section I VII 



II. DISTRIBUTION OF OCCUPIED STATES IN THE NEWNS-ANDERSON MODEL 

Our expression for the distribution of electronic states is obtained from the one-electron density matrix and is 
derived within the framework of the time-dependent, mean-field Newns-Anderso n 18 ' 19 model. This is defined by the 
Hamiltonian 

H(t) = ^H a {t) -Un aa {t)n a - a {t), (la) 

a 

with 

H a {t) = eaaWCcaa + ^^cL^ + I^^afeWcLc^+H.C.) (lb) 

k k 

where H{t) represents the total energy of the system and H a {t) is the many-electron Hamiltonian for spin a. c' a(T 
and ct are the creation operators for electrons in the adsorbate, \a a ), and metal, \k a ), states, respectively. V a k(t) is 
the interaction potential, U is the intra-adsorbate Coulomb repulsion energy and is the energy of the metal state 
| fccr ) . e a<7 is the mean- field energy level of the adsorbate state and is defined as 

=e a {t) + Un a ^(t), (2) 

where the time-dependent occupation of the \a a ) state 



n aa {t) = {a a \h llT {t)\a a ) 



(3) 



3 



is determined by the one-electron density matrix 

{b' a \h la (t)K) EE ((4(t)cva(*)». (4) 

Here \b a ) refer to one of the basis states \a a ) and \k a ), ((.}) denotes a thermal average and the annihilation and creation 
operators are given in the Heisenberg picture. As in our previous paper— we assume that the evolution of the system 
is driven by the variation in the bare adsorbate energy level, e a (f), and the interaction potential, V a k(t)- The Coulomb 
repulsion energy U is assumed to be constant. In order to model the behaviour of a real system, the variation of these 
parameters can be estimated from DFT calculations, as discussed in our previous work-ii. 

Since the mean-field Hamiltonian is quadratic in the annihilation and creation operators, the time-evolution of the 
one-electron density matrix, n\ a (t), can be obtained by considering the time-evolution of the one-electron states of 
the one-electron Hamiltonian, 

K{t) = e aa (t)\a a ){a a \+^e k „\k a }(k a \ + ^(V ak (t)\a a ){k„\+li.c.). ( 5 ) 

k k 

We will use \fj, a (t)) to represent the time-evolving set of electronic states of the one-electron Hamiltonian h a . These 
will evolve from an initial state at time to which is one of the basis states \b a ) of the system (i.e. \a a ) or \k a )). We 
note here that the initial interaction potential must be zero (i.e. V a k(to) — 0). The time dependence of the creation 
and annihilation operators is now simply determined by the time dependence of the one-electron states as 

M*)=£><t|mU*)>W*o) (6) 
v 

where \^' a {to)) — \b'). This result and (H]) shows that fi\ a {t) can be expressed in terms of \n a (t)) and the initial 
occupation /(e^ 4 ) of the dynamical state \[i a {t)) at time t$, where / is the Fermi function, as 

ftia(0=£M*)>M*)l/(&,)- (7) 

The quantity we are interested in is the distribution of occupied one-electronic states and how this evolves with 
time. In order to derive an expression for this distribution function, we need to take the diagonal matrix elements of 
h\ a {t) with respect to eigenstates of h a . We label these eigenstates as \v at ) with energies e£ t , and the distribution 
function becomes 

n a (e,t) ee ^(f<rt|rci<rKt)£(e-e£ t )> 

= J2\(v«tW(t))\ 2 fK t0 )S(t-<t)- (8) 

It is important to note that we choose states \v a t) that are not the usual eigenstates of the Newns- Anderson model. In 
the static case e a (t) and V a k (t) are held constant at a given value and (|lb|) and ([2]) are solved self-consistently, which 
means that n aa and e aa are determined entirely by these parameters. This is what we refer to as the adiabatic state 
of the system. In our system, however, n aa is not the self-consistent solution, but is the occupation of the adsorbate 
orbital at a given instant and therefore depends on how the system has evolved. The states \v a t) are the instantaneous 
states of the one-electron Hamiltonian rather than the adiabatic states. This choice of eigenstate means that the total 
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energy of the system is given by the first moment of the distribution function, that is 

E{t) = J2 I deen a (e,t) - Un aa (t)n a _ a (t) (9a) 
= Y.^AtW^WWi^fKJ-Un^na^it) (9b) 

<7,/X 

= YsdHAt))) -Un aa (t)n a _ a (t). (9c) 

cr 

This relationship between the distribution function and the total energy provides a useful check on the result for 
n a (e,t) presented later. We also note that in the long time limit, when evolution of the system is finished, the 
eigenstates \v at ) converge to the adiabatic states. 

In order to express the distribution function in a more manageable form we use Green's functions to replace the 
instantaneous states \v at ). We rewrite (jHJ as 

rvM) = --Im{Tr[n lCT (i)G ff (e;i)]}, (10) 
where G a is the instantaneous Green's function defined as 



\Vat){Vg 

-c + 

with r\ a positive infinitesimal. By introducing the basis set \b a ) into (|10p we find 

(12) 



VM) = -ilm|^n66'«r(*)G^,(e;t)| 



where G% b , = (b a \G a \b' a ) and nw a = (fo (J |ni (T (t)|6^}. n aa (t) = n aaa (t) is the dynamically evolving adsorbate occupa- 
tion which appears in (|la|) and ([2|) . In our previous work-i£ we obtained expressions for n aa (t) , but here we also need 
the occupation functions n aka (t), n kaa {i) and n kk i a (t), as well as the full set of instantaneous Green's functions. 
The Green's functions can be found from the Dyson equation 

G(e; t) = G°(e; t) + G°(e; t)V(t)G(e; t), (13) 

where G° is the unperturbed Green's function, and V is the interaction potential. For our system we have 

G1M) = CMl + CMlE^O^M), (14a) 

k 

G a ak (e;t) = G a :(e;t)J2Vak'(t)Gl, k (e;t), (14b) 
k' 

GLfe*) = G° k l(e;t)V: k (t)Gl a (e;t), (14c) 
G% k ,(e;t) = GlUe;t)5 k , k ,+GlUe-,t)V: k (t)G: k ,(e-,t), (14d) 

where 

Gll(e;t) = _ 1 , (15a) 

e - e aa [t) + irj 

G° k %(e;t) = (15b) 

e - e ka + if} 



5 



In order to find solutions to these equations we make two standard assumptions 1 -^^. First we assume that the inter- 
action potential V ak can be separated into a complex constant and a real, state-independent, time-varying function. 
Second we assume that the resonance width T, defined as 



r(t) = 27rJ2\Vak{t)\ 2 S{e~e k<J ) 



(16) 



is independent of energy e. This is often referred to as the wide-band limit. Using these approximations (|14ll becomes 



GL(e;t) = 
Gl k {e;t) = 
Gl a (e;t) = 



1 



e - e aa (t) ' 
1 



V ak (t) 



(e - e aa (t)) ' (e - e ka + irj)' 
(e-e ao .(f))'(e-e feo . + i7?)' 



where 



e - efeo- + ill ( e ~ £ fccr + irj) ' (e - e aC r(i)) ' (e - efc' CT + irj) ' 



(17a) 
(17b) 
(17c) 
(17d) 

(18) 



These Green's functions are very similar to those obtained by Anderson^, with the exception of the use of the 
instantaneous energy level e a o-(t) rather than the adiabatic level. 

Substituting (JTTJ) into (|12p . and noting that n aka — rit , allows us to write 

n a (e, t) = n a<T {t)p^ st \e, t) + ^ n kkrT (t)S(e - e ka ) 



-Im ■ 



1 



7T 
1 . 



-Im< 



e - e aa (t) 
1 



de' 



e - e c 



r{t) 



e — e + irj 
de 1 



.2Rc 



^ n aka (t)V ak (t)8(e' - e ktT ) 
L k 

de" 



e — e + it) J e — e + it) 

iJ2 n kk'At)V: k (t)V ak ,(t)S(e' - e ka )5{e" - e k , a ) 1, 



(19) 



where we have introduced energy integrals to make the recovery of the width T(t) easier, pia^ {e,t) is defined as 

r(t) 



27r[(e-e QCT (0) 2 +r(t) 2 /4] 



(20) 



which is very similar to the familiar adsorbate projected density of states, with the exception that once again the 
energy e aa (t) is the instantaneous rather than the adiabatic level. 

In order to complete our expression for n a we require n aa , n aka and n kk i a . Previously, we calculated n aa by solving 
the equation of motion for the creation and annihilation operators in the Heisenberg picture^. This solution is given 
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by (we assume H = 1) : 



dt\ exp 



+ exp 



i / e a „{t')dt' 



i / e aa (t')dt' 



^2 V a k{h) exp [-ie k(T (ti - i )] c kcr (t 



la 



from which we found 



* / d *i^afc(*i) ex P H e fer(* ~ *i)] c ao .(*i) + exp [-ie ktr (t - t )} c ka (t ), 



- f\{t')dt' + J def(e)\p a (e,t)\ 2 , 



n aa (t) = «4a(*)c a a(*)» = ^a CT (^o)exp 

where n aa (to) = f(^aa(to)) is the initial adsorbate occupation and p a is defined as 

Pa(e,t) -- 



(21) 
(22) 

(23) 



aii\/ — - — exp 



2tt 



i / (e a a(t') - e)dt' 



(24) 



As shown in Appendix [S] the remaining occupation functions can be evaluated in a similar manner using (|2ip and 
([22]) . Alternatively, the result for the occupation functions can be derived from ([7]) and solving for the time evolution 
of the states. 

The distribution function n a now follows from (|19p and the occupation functions in Appendix [X] The algebra is 
lengthy, but straightforward, with the definition of T(t) (JTSJ) used where necessary. We find 

n a (e,t)= [ de'f(e') q„(e, e', t) + p^ nst) (e, i)p CT (e', t) * - 2/( e )Re{g CT (e, e, i)} 



de'f(e')lm{p a (e\t)}Re 



r a (e,t)+p^ nst \e,t)exp 

' P [ : nst) (e,t) 



T(t')dt' 



to 



2 V^/ de ' /(e ' )Im ' 

J de'f(e')Re ■ 



1 m 
IT V 2tt 



e — e + ir\ 

q *(e,e',t)p ( r st) (e,t) ' 
e — e' + irj 

P { : nst) (e,t) \ 

(e-e' + z^) 2 f ' 



(25) 



where r CT and pa™ 11 ^ are defined as 



q a {e,e',t) = \ 

Jtn 



dh 



2tt 



P(T (e',ti) exp[i(e - e')(*i - t)], 



r a (e,t) = exp 



T(t')dt' 



la 



' , /r(*i) 

dt\\j — — exp 



2tt 



i I (e a(7 (t') - e)dt' 



ti 



r(t) 



27T e - e a(7 (i) ' 



(26) 
(27) 
(28) 
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such that, from ([20]) . 

(e,f) = |p^ st) (e,t)| 2 . (29) 

Each of these quantities, along with n aa and p a defined earlier, can be calculated numerically, as will be described in 
section HVl Given the variation of e a and T, and the value of U, the evolution of n a can therefore be computed. 

There are two tests which we can perform on the distribution n a to confirm that it is correct. First we can check 
that the number of electrons is conserved throughout the evolution of n a by integrating over all energies e; in appendix 
[5] we demonstrate that ([25]) obeys charge conservation. The second test is to calculate the total energy of the system 
by taking the first moment of the distribution function (|9a|) and comparing this to the rate of change of energy derived 
previously^. The verification that these approaches give the same result is long winded and will be presented in a 
future publication 20 . 



III. THE EXCITATION SPECTRUM 

In order to analyse the excitations of our system it is important to recall the definition of the distribution function 
n a (e,t); it is the time-evolving distribution of occupied electronic states. However, the quantity we require is the 
spectrum of excitations, for which we need to subtract an underlying distribution in which there are no electronic 
excitations. We write this instantaneous distribution as n CT t(e); it is given by ([8]) with \fx a (t)) = \v a t), which gives 



(e)=£/(C )5(e-C). (30) 



The difference between the times that appear in this expression is significant. f(e„ to ) is the initial occupation of the 
state that is connected to \v a t)', this occupation does not change with t. In a finite system, however, the eigenvalues 
can vary. In a system of N electrons, the change in the eigenvalue Ae^ t = — e^ to will be of order 1/N. We expand 
the Fermi function in (|30p to first order in Ae^ t , giving 



™<t<( £ ) = 



/(4 4 )-ac|(C) 



= /( e )^5( e -C)-|( e )^AC5(e-e^). (31) 

V V 

When integrated to give the total number of electrons in the sytem (as in Appendix B), the first term in this expression 
gives a quantity of order N electrons, while the second term integrates to a charge of order 1 electron (the sum over 
v contains N terms each of order l/N). In order to have an underlying distribution that conserves charge, it follows 
that this second term cannot be neglected in the limit as N — > oo. Higher order terms in the expansion of the Fermi 
function will yield of order 1/N electrons or less, and in the N — ► oo limit these terms can be ignored. 

The first term in <|3 1 1) can be dealt with using the Green's functions derived in the previous section. By introducing 
the basis states \b a ) the first term in (13T1) becomes 



m E k&*m i 2 ^ - o = -^E Im t)} ■ ( 32 ) 

t — ' 7T * — ' 

b,v b 
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When combined with (I17a|) and (|17d[) . and using the wide-band approximation, this yields 

f(e)p^(e, t)+Y / f(e k *)6(e - e ka ), (33) 
fc 

and therefore 



n,t(e) = f(e)pt st) (e, *) + £ - ^) - f £ A « e ~ ( 34 ) 

k 



A more convenient form for the third term in (fM)) can be determined by consideration of charge conservation. As in 
Appendix B, the integral of the instantaneous distribution function over energy: 

J den at (e) = J def(e)p^ (e,t) +£/(e te ) - / de| £ Ae^(e - C), (35) 

fc ' f 

should give the total number of electrons of spin cr in the system. The second term here gives the initial number of 
metal electrons, and so to conserve the number of electrons in n CTt (e) we must have 

-/ def e Y / Ae» t 5(e-e» t )=n a( ,(t )- J def(e)p^ st \e,t). (36) 

We now assume that the sum over the states \v G t) in (|3"6"|) is independent of energy over a range of several ksT either 
side of the Fermi level. This assumption is consistent with the wide-band approximation (see ((H])) used throughout 
this work. Equation (|36| then gives 

^Ae^(e-e^)=n OCT (t )- / def(e)p^(e,t) (37) 

and consequently the instantaneous distribution function becomes 



n„ t (e) = f(e)p^ st \e,t)+J2f(^)6(e-e kl7 )-f e n aa (t ) - J de'f(e')p^ (e\ t) 

k 

The difference between n CT (e,t), ([25)) , and n CT t(e), ([55)) , is the required spectrum of excitations nS e ^(e,t) 



(38) 



IV. COMPUTATIONAL METHODS 

In this section we outline the methods we have used to compute the distribution functions for the time-evolving, 
(I25p . and instantaneous, (|38[) . systems. To compute n a (e,t) we require three quantities; p ai q a and r a . As in reference 
liTJPcr, ([2"1)> . is obtained from 

^( e ',t) = -i{e a „{t) - 6'K( £ ',t) + (39) 

with the initial condition p<j(e', to) — 0. This equation is integrated on a finite grid of e' points using the fourth order 
Runge-Kutta method. The evolution of the system is driven by variation of e a and T with time, with the initial 
condition that T(t ) = 0. The energy e aa ([2]), and hence e aa (fl~8|) . is calculated using n aa which is obtained from (f23j) . 
A sufficiently fine e' grid is required to ensure the accuracy of the energy integral in (|23[) . q„ is obtained from p a by 
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using 



q*(e,e',t) = exp[-i(e - e')(t - t )] £ dh^^p^e' ,h) exp[i(e - e')(*i - to)]. (40) 

The e variable here is represented by a set of energy values at which the distribution function n a (e,t) is required. 
The time integral in (|40p is performed numerically using Simpson's method for each (e', e) grid point combination. 
The removal of any t dependence from the integrand in (|40|) makes it possible to evaluate q a at regular points in 
the Simpson's integration, rather than at a predefined time t. We use this property to explore the evolution of the 
distribution function n a , as e aa and T vary. r CT , (|27|) . is found in a similar manner to p a from 



— (e, t) = i(e aa (t) - e)r CT (e, t) + y exp 



r(t')*' 



to 



(41) 



with the initial condition r CT (e,to) = 0. The integration of r CT only needs to be performed for the energies e. 

Using (f39|) . ((40]) . (|4Tj) and the definition of n aa ([23]) . the first three terms in n a (e,t) can be calculated. The fourth 
term in n CT cancels with the second term in the instantaneous distribution function n at , (|38p . and can therefore be 
ignored. The final three terms in (|25p require further work due to the singularity in their integrands in the 77 — ^ + 
limit. 

The e' integral in each of the final three terms is separated into three sections; a window of width 2a around e' = e, 
a section below e' = e — a and a section above e' = e + a. The sum of the outer sections is similar to a principal 
value integral and we therefore use the notation PV Q to denote this. We then Taylor expand the integrand in the 
central window and perform this section of the integral analytically. We will demonstrate this 'analytic window' 
approximation using the fifth term in (|25| . ni , and then state the results for the sixth and seventh terms, and 

(7) 



The e' integral in n§ becomes 



™< 5) M) 



^Rc {^(e, t)} PV Q J de'I^llm {p a (e', t)} 



•- / de'f(e')lm{p a (e',t)}Re { P ° ^ , (42) 
7T J e _ a e - e' + n-j 



where we have dropped the rj from the first term as the PV Q integral does not cover the region in which it is important. 
We now expand the product f(e')p cr (e', t) as a Taylor series about e' = e to first order, yielding 



4 5) M) 



^Re (e, t)}pV a J de^lm {p a {e\ *)} 

+ -Re \p^ (e, t) d f } /(e)Im {p a (e, t)} 

+ -Re(^" s *)(e,t) f +a de' e ~/ ) ± \f(e)lm{ P<T (e,t)}] . (43) 
tt I. J e _ Q e-e' + «r?J del i 
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The integrals in the second and third terms in (I43|) can be evaluated and taking the — s- + limit we obtain 



+2 J-^(e e aa (t))p^(e,t)f(e)Im{p a (e,t)} 



-teyJ?S P £"He,t) (|lm{p <r ( C> t)}+/(e)Im{^(e ) t)}) , (44) 

where we have used the definitions of p%? at) , (El, and pi" lst) , IJ2SJ. The energy derivative of p a is obtained using a 
finite centred difference method. We note here that this expansion includes all first order terms in a, and will give 
the exact result in the a — > limit. 

By using the same method for we find 

nW(e,t) = 

+ (2tt - 4a{t - to)) V^/(^)Re {<£(e, 6, t)p^ (c *)} 

+4 °\/Wl im {' ?:(e ' e '^ mst)(e ' t) } 

+4a \/W /(e)im |^ mst) ( £ ' *) jf (¥ (e ' ti} + z(ti " *° K(e ' <} ) } • (45) 

The final term in (|25p . ni 7 \ requires a little more work before applying our approximation due to the l/(e — e') 2 



dependence of the integrand. Integrating ria by parts yields 



-/(e') ] £ + , y 4f 1 



e — e' + z?/ 



de' e — e' + irj 



_ _r(0.f"M - [A 1 , (46, 

2vr (e-e'_) tt V 2tt J de' 1 e - e' + irj J v ; 

where e' + and e'_ are the upper and lower limits of the e' range over which we are integrating, and we have assumed 
that f{e'i) = and f(e'_) = 1. The remaining integral in (|46p can be dealt with in the same manner as and , 
giving 

2ir (e — e_J ztt J dee — e 

-(e - e a 4t))f/-He,t) + £^^p(^») (e, t). (47) 

This result, like ([4"4")l and (|43|) . includes all first order terms in the window half- width a. 

The instantaneous distribution function is more straightforward to calculate than n a . The first term in n a t, l|38[). 
can be calculated directly using the definition of p£a St \ (f20|) . The second term in (|38|) cancels with an identical 
term in n CT , (I25|) . The final term in (|38[) cannot be used directly due to the truncation of the e' range over which we 
perform numerical integrals. We require the integral of n a t over all energies e to conserve the number of electrons and 
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de 



consequently we modify the final term in (|38[) to 

/>oo 

Tw(to)-/ &/(e)p£ lrt )( e ,t) 

J — oo 

u OCT (t )- fV/(e')p£r'V') 

na<r(io) 



4f 



' de 



(48) 



where we have assumed f(e' + ) = as before and have performed the integral over e' up to e'_ analytically. 

It is convenient when performing numerical calculations to work with de-dimensionalised parameters. As in our 
previous worki£ we scale all parameters by a width To; so that, for example, tjj = Tot, = T/ro, n a D — ron CT , 
= e/r , f/o = U/Yq and fc^T^ = kgT/Y where we have used subscript -D to denote the dimcnsionlcss quantities. 
Typically in practice To will have a value between 1 and 3 eV. 

We have carried out extensive tests of the numerical stability of our calculations. We find that for run times of 
up to 200 dimensionless time units, well converged results are obtained for 160,001 e' D points in the range —20 to 
20, where the Fermi level is fixed at e' D = 0. The analytic window half-width a is set to be the same as the e' D grid 
spacing and the Runge-Kutta integration step is 0.01 dimensionless time units. The use of a non-zero temperature 
in our model aids the stability of n a D, particularly in the region around the Fermi level. In this work we have used 
a value of UbTd = 0.02, which typically corresponds to a temperature of the order of a few hundred Kelvin. With 
these parameters we have found that our computed distribution conserves charge over the full length of a dynamical 
run to better than 10 of an electron. 



V. NUMERICAL RESULTS 

We have chosen three model systems to demonstrate the non-adiabatic behaviour of the Newns- Anderson system. 
In each case Ud = 3 and e a D is held fixed. The model is driven by varying Yd from to 3 over 50 dimensionless 
time units using an error function with a peak gradient of dYo/dtD = 0.3 at tn = 25. The three model systems 
we have chosen differ in their bare energy levels t a D\ we use e a D = —2.5, —1.5 and —0.5. Each of these systems is 
driven through the spin transition by the Yd variation. The first system has the minority spin energy level closest 
to the Fermi level ep , the second has the energy levels e OCT equidistant from tp and the final system has the majority 
spin energy level closest to cf- These parameters were chosen to be broadly comparable with those extracted from 
DFT calculations of the H/Cu (111) system explored previously^, with the exception that the rate of change of the 
parameter Y has been increased to exaggerate the non-adiabatic behaviour. 

Figure [1] shows the evolution the adsorbate state occupations and the mean-field energy levels, along with the 
adiabatic equivalents calculated as in ref 17, for each of these model systems. The e a D = —2.5 calculations show 
the minority spin energy level crossing the Fermi level, gaining occupation and converging to the majority level. 
The e a D — —1.5 levels converge simultaneously on the final state, with small deviations from the analytic result 
{n a <j = 1/2, e aa — ep) due to the truncation of the e' D range in the numerical calculations. The final model system, 
with e a D — —0.5, involves the majority level crossing the Fermi level and approaching the falling minority level, 
resulting in a low occupancy final state. In each of the model systems the adiabatic occupations exhibit a sharp 
transition from a spin-polarised to a spin-degenerate state at around to = 25. The dynamical occupations, however, 
overshoot this sharp spin transition and the net polarisation then falls roughly exponentially to below 0.01 by to — 35. 

Figures [3] and |4] show a series of snapshots of the evolution of the electronic excitation spectra for each of the 
model systems. In each case the majority of the evolution occurs during the period in which the rate of change of Yd 
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FIG. 1: Adsorbate state occupation and energy level variation for the adiabatic (dashed lines) and time-evolving (solid lines) 
models: (a), (c) and (e) occupations for e a D = —2.5, —1.5 and —0.5 respectively, (b), (d) and (f) mean-field energy levels for 
e a D = —2.5, —1.5 and —0.5 respectively. The dotted line in (b), (d) and (f) denotes the Fermi level. The arrows in (a) and (b) 
denote the spin-states, with f indicating majority spin and J, minority spin. 

is largest. The snapshots of the e a D = —2.5 system, in Fig. [2j show the early evolution of occurring primarily in 
the minority spin-state, while the majority state evolves later. The e a D = —1.5 model system in contrast, see Fig. [3l 
has both spin states evolving in a symmetrical manner resulting in near-identical spectra for majority spin electrons 
and minority spin holes. The e a u = —0.5 model system, see Fig. [4] is similar to the e a o = —2.5 system with the 
reversal of the spin states; evolution occurs early and late in the majority and minority states respectively. 

In each of the model systems the excitation spectra show a number of similar features, with the balance between 
them determined by the parameters of the system. The majority spin spectrum consists of a fairly narrow peak of holes 
close to the Fermi level with a relatively small tail extending below ep . Above the Fermi level, however, the majority 
spin distribution is flatter and broader, leading to a larger tail of electronic excitations at higher energies. The overall 
size of the majority spin distribution is governed by the total change in the occupation of the relevant adsorbate state. 
The largest change (from 1.00 to 0.36) occurs for e a D = —0.5 (Fig. [IJe)) and this gives rise to the largest majority-spin 
excitation spectrum (Fig. |4|). Conversely, the small change in the majority-spin occupation for e a D — —2.5 (from 1.00 
to 0.61, Fig. [UJa)) yields a relatively small excitation spectrum (Fig. These trends are reversed for the minority- 
spin components. These have a narrow electron peak near the Fermi level with a small number of higher-energy 
electron excitations. There is a broader hole distribution, which gives the dominant contribution to higher-energy 
holes below the Fermi energy. The magnitude of the minority-spin spectrum is again governed by the total change in 
the relevant occupation, so that the largest minority spectrum now occurs for e a o = —2.5. 

These results have interesting consequences in the context of chemicurrent generation, where electrons or holes with 
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FIG. 2: Excitation spectra snapshots for e aD = -2.5 at times (a) to = 22.00, (b) to = 23.25, (c) to = 24.50, (d) t D = 25.75, 
(e) t D = 27.00, (f) t u = 28.25, (g) t D = 29.50, and (h) after F D variation has finished (to = 50.00). Solid lines denote the total 
electronic excitation spectrum, dashed the spin T (majority) component and dot-dashed the spin J, (minority) component. 



sufficient energy to cross a Schottky barrier are detected 5 -. Figures [3] to H] show that an asymmetry in the adsorbate 
energy levels with respect to the Fermi level will lead to an asymmetry in the measured electron and hole currents. 
For example, the ratio of electrons to holes at cd = ±0.5 is 1:6.6, 1:1.2 and 5.0:1 for e a D = —2.5, —1.5 and —0.5 
respectively. We also note that the high energy tails of til^v are dominated by the majority spin for electrons and the 
minority spin for holes. For each of the model systems the spectrum consists of at least 96% majority spin electrons 
above = 0.2, with similar fraction of minority spin holes below — —0.2. This implies that a spin-polarised beam 
of adsorbates made incident on a metal surface would generate a spin-polarised chemicurrent. 

Figure [5] shows the effect of changing the rate of variation of Tjj on n£p . The rate of variation can be interpreted 
as the approach speed of an adsorbate, where larger peak gradients of imply higher speeds. In Fig. [5^a) we have 
plotted the total electron excitation spectrum for e a D = — 1-5 for three different speeds with a logarithmic scale for 
n ^D ■ F° r eacn °f the approach speeds the excitation spectrum above to ~ 0.2 varies exponentially with energy. As 
would be expected, decreasing the speed of approach to the surface reduces the magnitude of i.e. the evolving 

distribution of occupied states will converge eventually to the instantaneous distribution. To analyse this further, we 
have fitted the electronic excitation spectra above €d = 0.2 to the exponential e XtD for a number of different approach 
speeds. In Fig. E^b) the decay parameter A is plotted as a function of approach speed for the three model systems. We 
find that the variation of the parameter A is well modelled by a speed -0 ' 5 dependence for each of the model systems, 
with the E a £> = —1.5 system having a larger gradient than the e a u — —2.5 and —0.5 systems, which behave similarly. 
We have not, to date, been able to explain the origin of this speed -0 5 dependence. 
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VI. CONCLUSIONS 

We have presented an analytical solution for the time evoution of the mean-field Newns-Anderson model, which 
has enabled us to calculate the spectrum of hot electrons which are excited in the course of an encounter between 
an absorbate and a metallic surface. Although the Newns-Anderson model is a grossly over-simplified description of 
any real system, it does have the major advantage that results can be obtained quickly and that trends can be easily 
investigated. In this paper we have focussed on model systems that have a spin transition; in a future publication 
we will examine in detail the excitation spectra for H/Cu, H/Ag and O/Ag, all of which exhibit this transition 21 . 
However, our analysis can be applied to a wide range of other systems; all that is required is a model for the variation of 
the bare energy level and the resonance width. The Newns-Anderson model should provide a useful semi-quantitative 
description of the non-adiabatic coupling between an absorbate and a metallic substrate. This will allow, for example, 
an investigation of the validity of a nearly adiabatic (ie friction-based) approach and the extent to which the forced- 
oscillator model provides an accurate description of the excitation spectrum. The Newns-Anderson model will also 
provide a useful check on a more sophisticated theory of electronic excitations, for example, one based on ab-initio 
time-dependent density functional theory 22 . 

There is one aspect of our work which merits further discussion. One of the key advantages of using a simple 
analytical model is that the results obtained from it are usually transparent, which enables us to obtain physical 
insights that can be applied to more complex situations. In our case, however, the expressions for the time evolution 
of the electronic distributions become so complicated that this transparency is lost. Although it is possible to determine 
where each term in the final expressions have come from, we have not to date been able to extract a clear physical 
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FIG. 5: (a) Electron excitation spectrum at to = 50 for e a D = —1.5 at full speed (solid line), half speed (dashed) and quarter 
speed (dot-dashed line) relative to the calculation in Fig. [3] (b) variation of the decay parameter A (see text) with speed for 
e a D = —2.5 (dashed line and diamonds), e a D = —1.5 (solid line and circles) and e a D = —0.5 (dot-dashed line and squares). 
Points are fitted gradients to n^D an d h nes are speed -0 ' 5 fits to these data. 



16 



picture of the excitation process. Clear trends can be observed in Figs (2) to (4), and in Fig (5), but going from a 
description to an explanation remains a challenge. 
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APPENDIX A: OCCUPATION FUNCTIONS 



Substitution of (21]) into (22|) yields 



t rti 

diiV r a * fe (ii)exp[-ie fe(T (t-ti)] / eft 2 exp 

to J to 



-i \ e aa {t')dt' 
Jt 2 



-i / dtiV* k {ti)ex.-p[-ie ka (t -ti)}exp 

Jto 

+ ex.p[-ie ktT (t - t ))c ka (t ). 



^ V ak > (t 2 ) exp [-iek'a(h ~ to)] e feV (*o) 

C a a(*o) 



-i / i aa {t')dt' 
Jt 



(A.l) 



The required occupation functions follow by using (|2"Tj) and (|A.1|) together with the initial conditions 

((dUhKM)) = f(e a M), ((cL(to)£ k M)) = 0, ((iM&vM)) = fM6 kk >. We hnd 



n aka {t)= ((4rWc fcCT (i)» 



-i dhVMh) I de'f(e')p*(e',t) P<T (e',t 1 )eMi{e'-eka)(t-ti)] 

Jto 



+if{e k a) / dtiV7 fc (ti)exp 
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* / (eL(t') ~ £k<r)dt' 
ti 



(t ) exp 



x / dtiV* k (h)exp -i ( (e a<T {t') - e ka )dt' 



(A.2) 



and 



n kk , a (t) = ((c{jt)c k .(t))) 



= / de'f(e') ^dt 1 V r ait (ti)p*(e',t 1 )ex P [-i( efcff - e ')(ti-t)] 
x / dt 2 V;V(i 2 )p (7 (e' ) t 2 )exp[i(e^ ff - e')(t 2 - t)} 
-/(efev)/ dtiVofc^i) exp[i(e fc(T - e^X* - 
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/ *lKfe'(*i) exp[i(e fc(T - e fc , (T )(i - h)] 
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to 



x / dt 2 V* k ,(t 2 )exp 

Jto 



-i / (eaa(t') - e k/a )dt' 

Jt 



f(tkcr)5kk' 



The remaining occupation function in (|19|) is n kka , which is a special case of n kk > a . From (IA.3|) we find 
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(A.3) 



n kklJ {t) = / def(e') 



dtiV£ k (ti)p (T (e',ti)exp[i(eit<T - e')(*i _ *)] 



to 

-2/(e fcCT )Re| / dhV a * k (h) [ ' dt 2 V ak (t 2 )exp 
I -/to -'to 



Jt 2 



+«aa(*o) 
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(A.4) 



APPENDIX B: CHARGE CONSERVATION OF THE DISTRIBUTION FUNCTION n a 



In this appendix we demonstrate that the time-dependent distribution function (|25p conserves charge by taking the 
integral of n a over all energies. We will use numerical superscripts to denote the individual terms on the right hand 
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side of (|25|) . The integral over the first term can be written, using (f26|) and {28)), as 



den«(e,t) = de 1 f{e') de q a (e, e',t) + p% nst \e,t)p a (e' ,t) 



= lde'f(e') dh / dt 2 vT (t i ) r (i 2 )p* (e' , < 2 K (e' , t x ) exp [id (t 2 - 1 1 )] 

'to "'to 

exp[ie(ti - i 2 )] 
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x de 



e~xp[—ie(t — ti)] 



+ / de7(e')b ff (c',t)| / * 



r(i) 



'2^[(e- £ - QCT (i)) 2 +r(t) 2 /4]' 



(B.l) 



The e integral in the first term of this expression gives the delta function 6(ti — i 2 ) and that in the third is unity. 
The e integral in the second term can be evaluated using contour methods, with the contour closed in the lower half 
plane. By the residue theorem this integral is zero, because the pole at e = e* aa {i) — e aa {t) + iT{t)/2 is in the upper 
half plane. Equation (|B.1[) therefore simplifies to 



den^(e,t) = jf faFfo) J de'f{e') \p„{e' M)\ + J de' f(e') \ Pa (e' ,t)\ 2 . (B.2) 
The integral of the second term in (|25| results in a term which we cannot simplify at this point so we simply state 

(B.3) 



denf\e,t) =-2 1 dh 
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The third term in (|25|) can be integrated to give 



den^\e,t) = n aa (t ) / de 



r rT (e,t)+p^ nst He,t)exp 
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Using contour methods, and the definition of r a (I27p . this expression can be evaluated, yielding 



denf>{e,t) =n aa (t ) / faTfe) 



1 exp 



T{t')dt' 



n QCT (t )exp 



T(t')dt' 



The fourth term in (|25| will give us the number of metal electrons, 7V™ eta ': 

J den^(e,t) = J de^ fMS(e - e ktr ) = £ /( 



e k a) = N, 



t 



metal 



(B.5) 



(B.6) 



The fifth, sixth and seventh terms in (|25p can easily be shown to integrate to zero using contour methods. 
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By combining (|B.2[) . ()B.3|) . (|B.5p and ()B.6|) . and using the definition of n aa (t) in (|23|) . we obtain the following 
expression for the total charge; 



de n a (e, t) 
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+N metal_ ( R7 ) 

To remove the time integrals from (|B.7|) we consider the time derivative of n aa (t). Using (f23|) this becomes 



d , . i ^ d 

-r;n a(7 {t) = n aiy (t )— exp 
dt dt 



T{t')dt' 



+ J def(e)j t \p a (e,t)\ 
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"By integrating over time from £q to t this yields 



(t) - n aa (t ) 



dt 1 T{t l )n aa {t 1 ) + 2 / dh 



to 



which, on combination with (|B.7[) . gives 



(B. 



J t *iy^ / def(e)Re{ Pa (e,t 1 )}, (B.9) 



de n<T {e,t) =n aa (t )+N2 



■metal 



(B.10) 



This confirms that charge is conserved, i.e. the number of electrons of spin a is the sum of those in the adsorbate and 
the metal prior to any interaction and does not vary with time. 
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